function ParetoOptimalTracing(nelx,nely,desVolFrac,problem)
% Ref: "A 199-line Matlab Code ...", K. Suresh, SMO, Volume 42, Issue 5, pp 665-679, 2010
if (nargin == 0)  % default values
    nelx = 40;nely = 20; % The grid size for topology optimization
    desVolFrac = 0.5;  % The final volume fraction desired
    problem = 1; % 1 or 2 for cantilevered beam problems
end
paretoAggressiveness = 0.65;%higher (upto 1), closer to Pareto curve, but may terminate early
volDecrMax = 0.05; % step-size for pareto-tracing
totalIter = 0;elemsIn = ones(nely,nelx); % intialize the domain 
volIndex = 1;volFractions(1:500) = 1; J(1:500) = 1;volfrac = 1; JTemp = 0;
[U,J(1)] = FE(nelx,nely,elemsIn,problem); % Solve FEA
T = ComputeTSField(U,elemsIn); % Compute topological sensitivity
while (volfrac > desVolFrac)
    volfrac = max(desVolFrac,volfrac-volDecrMax); % move to the next volume fraction
    localIter = 0;  success = 0;
    while (localIter < 10) % typically 2 ~3 is sufficient, but worst case     
        isParetoOptimal = isTopologyParetoOptimal(T,elemsIn,paretoAggressiveness);
        if (JTemp > 10*J(1)),break;end;% divergence
        if ((localIter > 1) &&(isParetoOptimal)), success = 1;break; end % success
        value = findCutoffWithVolumeFraction(T,volfrac);
        elemsIn = ones(nely,nelx); % start with the full domain
        elemsIn(ind2sub(size(T),find(T < value))) = 0; % remove elements
        [U,JTemp] = FE(nelx,nely,elemsIn,problem); % FEA
        T = ComputeTSField(U,elemsIn); % Topological Sensitivity
        localIter= localIter+1;totalIter= totalIter+1; 
    end
    if (~success),break;end;% done
    volIndex = volIndex+1; volFractions(volIndex) = volfrac; J(volIndex) = JTemp;
    value = findCutoffWithVolumeFraction(T,volfrac); % as above
    figure(1);plotContour(T,value,figure(1));
    title(['v=' num2str(volfrac) '; J = ' num2str(J(volIndex)) '; #FEA = ' num2str(totalIter)]); 
end
figure(2); plot(volFractions(1:volIndex),J(1:volIndex),volFractions(1:volIndex),J(1:volIndex),'*');
xlabel('Volume'); ylabel('Compliance'); grid on;
function [U,J]=FE(nelx,nely,elemsIn,problem)
if (problem == 1) % Cantilevered beam; 
    fixeddofs = 1:2*(nely+1); % left edge
    forcedDof = 2*(nelx+1)*(nely+1)-nely; % y force
elseif (problem == 2) % MBB beam
    fixeddofs = 1:2*(nely+1); % left edge
    forcedDof = 2*(nelx+1)*(nely+1); % y force
end
K = zeros(2*(nelx+1)*(nely+1), 2*(nelx+1)*(nely+1));% start with full matrix
F = sparse(2*(nely+1)*(nelx+1),1); U = zeros(2*(nely+1)*(nelx+1),1);
[KE] = elementStiffness();
for elx = 1:nelx
  for ely = 1:nely
    n1 = (nely+1)*(elx-1)+ely; 
    n2 = (nely+1)* elx   +ely;
    edof = [2*n1-1 2*n1 2*n2-1 2*n2 2*n2+1 2*n2+2 2*n1+1 2*n1+2]';
    K(edof,edof) = K(edof,edof) + max(1e-8,elemsIn(ely,elx))*KE;%to avoid singularity
  end
end
K = sparse(K); F(forcedDof,1) = -1;
freedofs    = setdiff(1:2*(nely+1)*(nelx+1),fixeddofs);
U(freedofs,:) = K(freedofs,freedofs) \ F(freedofs,:); U(fixeddofs,:)= 0;
J = F'*U;
function [KE]=elementStiffness() % element stiffness
E = 1.; nu = 0.3;
k=[ 1/2-nu/6   1/8+nu/8 -1/4-nu/12 -1/8+3*nu/8 -1/4+nu/12 -1/8-nu/8  nu/6  1/8-3*nu/8];
KE = E/(1-nu^2)*[ k(1) k(2) k(3) k(4) k(5) k(6) k(7) k(8)
                  k(2) k(1) k(8) k(7) k(6) k(5) k(4) k(3)
                  k(3) k(8) k(1) k(6) k(7) k(4) k(5) k(2)
                  k(4) k(7) k(6) k(1) k(8) k(3) k(2) k(5)
                  k(5) k(6) k(7) k(8) k(1) k(2) k(3) k(4)
                  k(6) k(5) k(4) k(3) k(2) k(1) k(8) k(7)
                  k(7) k(4) k(5) k(2) k(3) k(8) k(1) k(6)
                  k(8) k(3) k(2) k(5) k(4) k(7) k(6) k(1)];           
function [T] = ComputeTSField(U,elemsIn)
[nely,nelx] = size(elemsIn);T(1:nely,1:nelx) = 0; % all zeros
gradN =0.5*[-1 1 1 -1;-1 -1 1 1]; % at center
E0 = 1;nu = 0.3;D0 = 1/(1-nu^2)*[1 nu 0; nu 1 0;0 0 (1-nu)/2]; % plane stress
for elx = 1:nelx
  for ely = 1:nely
    if (~elemsIn(ely,elx)), continue; end; % skip outside elements
    n1 = (nely+1)*(elx-1)+ely;  n2 = (nely+1)* elx   +ely;
    edof = [2*n1-1 2*n1 2*n2-1 2*n2 2*n2+1 2*n2+2 2*n1+1 2*n1+2]';
    uGrad = gradN*U(edof(1:2:end));vGrad = gradN*U(edof(2:2:end));
    strains = [uGrad(1); vGrad(2); (uGrad(2)+vGrad(1)) ];
    stresses = D0*E0*strains;
    stressTensor = [stresses(1) stresses(3); stresses(3) stresses(2)];
    strainTensor = [strains(1) strains(3)/2; strains(3)/2 strains(2)]; 
    T(ely,elx) = 4/(1+nu)*sum(sum(stressTensor.*strainTensor))-(1-3*nu)/(1-nu^2)*trace(stressTensor)*trace(strainTensor);
  end
end
T = filter2(fspecial('gaussian',[3,3],0.6),T); % Smoothen the field
function value = findCutoffWithVolumeFraction(field,volfrac)
nElemsToRemove = ceil(numel(field)*(1-volfrac));
sortedField = sort(field(:));
value = mean(sortedField(nElemsToRemove));
function [isParetoOptimal] = isTopologyParetoOptimal(T,elemsIn,paretoAggressiveness)
T_InMin = min(T(elemsIn==1)); % Min of topological field inside the domain
T_OutMax = max(T(elemsIn==0)); % Max of topological field outside the domain
if (T_InMin > paretoAggressiveness*T_OutMax), isParetoOptimal = 1;else isParetoOptimal = 0; end
function plotContour(T,value,fig)
[nely,nelx] = size(T);
figure(fig);clf;fill([1 nelx nelx 1],[1 1 nely nely],'b'); hold on;
contourf(-T(end:-1:1,:),[-value -value]); axis('equal'); axis tight;axis off;